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We present a numerical code designed to conduct a likelihood analysis for clusters of nucleons above 
10" eV originating from discrete astrophysical sources such as powerful radio galaxies, 7— ray bursts 
or topological defects. The code simulates the propagation of nucleons in a large-scale magnetic field 
and constructs the likelihood of a given observed event cluster as a function of the average time delay 
due to deflection in the magnetic field, the source activity time scale, the total fiuence of the source, 
and the power law index of the particle injection spectrum. Other parameters such as the coherence 
length and the power spectrum of the magnetic field are also considered. We apply it to the three 
pairs of events above 4 x 10^^ eV recently reported by the Akeno Giant Air Shower Array (AGASA) 
experiment, assuming that these pairs were caused by nucleon primaries which originated from a 
common source. Although current data are too sparse to fully constrain each of the parameters 
considered, and/or to discriminate models of the origin of ultra-high energy cosmic rays, several 
tendencies are indicated. If the clustering suggested by AGASA is real, next generation experiments 
with their increased exposure should detect more than ~ 10 particles per source over a few years 
and our method will put strong constraints on both the large-scale magnetic field parameters and 
the nature of these sources. 
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I. INTRODUCTION 

Despite more than 30 years of experiments on ultra- 
high energy cosmic rays (UHECRs), i.e., air showers with 
primary energy E > 1 EeV {— lO^^EeV), the origin of 
these particles remains unknown. The most distinctive 
feature of the energy spectrum of UHECRs is undoubt- 
edly the "ankle", at £; ~ lO^^-^ eV The cosmic ray 
spectrum is steep up to the ankle, corresponding to a dif- 
ferential power law index of ~ —3.1, and flattens beyond, 
with an index ~ —2.7. Data from the Fly's Eye experi- 
ment also suggest that the chemical composition is dom- 
inated by heavy nuclei, presumably iron, up to the ankle, 
and by protons beyond Hence, the ankle might mark 
a transition from galactic to extra-galactic origin [|l],D. 
This idea is supported by the observed isotropy of cosmic 
rays for i5 > 4 x 10^^ eV [||, together with the fact that 
charged particles with E > 10^^ eV have too large a gyro- 
radius in the galactic magnetic field to be confined. How- 
ever, in analogy to the 7— ray burst (GRB) phenomenon, 
an origin of UHECRs in an extended galactic halo cannot 
be excluded thus far. 

If UHECRs with E > lO^^-^eV are indeed protons 
of extra-galactic origin, one would expect to detect the 
so-called Greisen-Zatsepin-Kuzmin cutoff |Q] (hereafter 
GZK) due to photopion production on the cosmic mi- 
crowave background (CMB) by nucleons of energy E > 



70 EeV. The presence of this cutoff has been suggested 
by various experiments, although the energy spectrum 
around E ~ 100 EeV is yet to be agreed upon. Nonethe- 
less, the combined data of the Haverah Park, Yakutsk, 
Fly's Eye, and AGASA experiments give 7 detected 
events above 100 EeV, for 22 expected from extrapola- 
tion of lower energy data, which provides a strong hint 
toward the presence of a cutoff j|. Finally, a 2cr de- 
tection of this GZK cutoff has been recently claimed by 
the AGASA experiment from the combined data 1991- 
1996 §. 

Photopion production on the CMB limits the range of 
nucleons above 100 EeV to about 30 Mpc, whereas heavy 
nuclei are photodisintegrated on an even shorter distance 
scale 1^. In this frame, the detection of UHECRs with 
E > 100 EeV (7[§,|,|,| has triggered considerable discus- 
sion in the literature on the nature and origin of these 
particles ]lO|-|l^. Indeed, even the most powerful astro- 
physical objects such as radio galaxies, quasars, and ac- 
tive galactic nuclei are barely able to accelerate charged 
particles to such energies [T^ ]. Moreover, there is no obvi- 
ous optical or radio counterpart to these UHEGR events 
at a distance less than a few hundred Mpc. For the high- 
est energy event observed, ~ 3 x 10^" eV [Q , the quasar 
3C 147 is the best candidate as far as energetics is con- 
cerned; however, it lies ~ IGpc away, and if the pri- 
mary were a proton, its injection energy would have to be 
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E > 10^^ eV. A similar problem arises for the less likely 
option of a 7— ray primary |^2|, whereas neutrino pri- 
maries in general imply too large a flux because of their 
small interaction probability in the atmosphere 

Currently, there are three main classes of models for 
the origin of UHECRs. The most conventional one as- 
sumes first order Fermi acceleration of protons at astro- 
physical magnetized shocks (see, e.g., Ref. ||l5[). This 
mechanism is usually associated with prominent astro- 
physical objects such Fanaroff-Riley Class II radio galax- 
ies, either in the hot spots or in the lobes |16| and could 
accelerate protons up to ~ 10^^ eV. 

It has also been suggested that UHECRs could be as- 
sociated with cosmological 7-ray bursts (GRBs) p7|-pO|. 
This is mainly motivated by the fact that the required 
average rate of energy release in 7— rays is compara- 
ble to the one in UHECRs above lOEeV turn out to 
be comparable. Protons could be accelerated beyond 
100 EeV within the relativistic shocks associated with 
fireball models of cosmological GRBs Since the 

rate of cosmological GRBs within the field of view of 
the cosmic ray experiments which detected events above 
100 EeV is about 1 per 50 yr, a dispersion in UHECR 
arrival times of at least 50 yr is necessary to reconcile 
the observed UHECR and GRB rates. Such a disper- 
sion could be caused by the time delay of protons ac- 
quired through their deflection in large-scale magnetic 
fields (LSMF) |p^ , ^ . However, this requirement could 
be circumvented in a model where UHECRs are produced 
in a galactic halo population of GRBs . 

Finally, so-called "top-down" models constitute an- 
other class of scenario. There, particles are created at ex- 
tremely high energy in the first place by the decay of some 
supermassive elementary "X" particle associated with 
new physics near the grand unification scale Such 
theories predict phase transitions in the early universe 
that are expected to create topological defects such as 
cosmic strings, domain walls or magnetic monopoles. Al- 
though such defects are topologically stable and would be 
present today, they could release X particles due to phys- 
ical processes such as collapse or annihilation. Among 
the decay products of the X particle are jets of 10^ — 10^ 
hadrons most of which are in the form of pions that subse- 
quently decay into 7-rays, electrons, and neutrinos. Only 
a few percent of the hadrons are expected to be nucle- 
ons ||2^. Thus, typical features of these scenarios are the 
predominant release of 7-rays and neutrinos, and spectra 
that are considerably harder than in the case of shock 
acceleration. For more details about these models, see, 
e.g., Ref. H. 

Recently, a possible correlation of a subset of events 
above 40 EeV among each other and with the supergalac- 
tic plane was reported by the AGASA experiment 
Among 20 events with energy above 50 EeV, two pairs of 
events with an angular separation of less than 2.5° were 
observed within 10° of the supergalactic plane at an an- 
gular resolution of ~ 1.6°. The probability for that to 
happen by chance for an isotropic, unclustered distribu- 



tion is ~ 4 X 10~^. A third pair was observed among 
36 showers above 40 EeV, with a chance probability of 
about 6 X 10^'^. Although the muon content of the show- 
ers suggests that the primaries are protons in each case, 
7-rays cannot be excluded in the following, however, 
we will assume that the primaries are protons. 

This suggests that the events within one pair have been 
emitted by a single discrete source possibly associated 
with the large-scale structure of the galaxies. The deflec- 
tion of a charged particle in a magnetic field is inversely 
proportional to its energy E. Therefore, the fact that the 
lower energy event in the pair with the greatest energy 
difference (see Table |) arrived later might hint to the 
pair's origin in a burst, i.e., on a time scale ^ lyr. The 
time delay would then be dominated by magnetic deflec- 
tion of the lower energy particle. Even the two other 
pairs observed by AGASA, in which the higher energy 
particle arrived later, could have originated in a burst. 
Since the average time delay due to magnetic deflection, 
te, scales as te oc E~'^, the dispersion in time delay 
is bounded from below by twice the energy resolution, 
At/t ^0.6, and could be as large as ~ 2 Further- 
more, the distance to the source cannot be much larger 
than ~ 100 Mpc, if the higher energy primary was ei- 
ther a nucleon, a nucleus, or a 7-ray, since its energy was 
observed to be ^ 75 EeV in all three pairs. 



Pair # 


Energy [EeV] 


date 






1 


210 


93/12/03 


0.4 


-41.1 




51 


95/10/29 


1.1 


-42.3 


2 


78 


95/01/26 


0.5 


56.9 




55 


92/08/01 


2.0 


55.3 


3 


110 


94/07/06 


57.3 


18.6 




43 


91/04/20 


57.8 


21.2 



TABLE I. Properties of the pairs observed by AGASA. 
6^*^ and 6^ are the supergalactic and the galactic latitude, 
respectively. 

Possible explanations of these observations rest not 
only on the model of UHECR origin but also on the 
largely unknown features of the LSMF (for a review see, 
e.g., ReL [^). Here, we call LSMF both extra-galactic 
fields with strengths lO'^^ G < B,ms < 10~^ G as well 
as fields in the halo of our Galaxy with 10^* G < Bh ^ 
10^^ G. In the present paper, we show how the distribu- 
tion of arrival times and energies of clusters of UHECR 
events originating from the same astrophysical source 
such as the ones suggested by AGASA can be used to 
constrain both the nature of the source and the struc- 
ture of the LSMF. We restrict ourselves to nucleons and 
perform Monte Carlo simulations of their propagation 
through the LSMF. Then we construct the likelihood as 
a function of several parameters characterizing the LSMF 
and the source of the pairs observed by AGASA. We de- 
scribe our Monte-Carlo simulations and the likelihood 
analysis in Sect. 2. In Sect. 3, we discuss our findings 
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and their implications for the structure of the LSMF and 
for the models of UHECR origin. Sect. 4 concludes with 
a discussion of the future prospects of the method pre- 
sented in this paper. 



II. NUMERICAL SIMULATIONS 
A. The Monte-Carlo code 

Previous works have focused on various aspects of non- 
diffusive, nearly straight-line propagation of extremely 
high energy nucleons in LSMF. In Ref. |^ the diffu- 
sion problem was studied, in the limits DOe ^ and 
DOe ^ ^c, where D is the distance to the source, is 
the r.m.s. deflection angle, which is a function of energy 
E, and Ic is the coherence length of the magnetic field. 
However, these authors did not consider pion production, 
and their results are only valid for energies E EeV. 
In Refs. [p8|-^ propagation codes that follow multiple 
species were constructed, taking all types of interactions 
into account, although neglecting any scattering on the 
magnetic field; thus these calculations cannot treat the 
time delay of UHECRs. Finally, the authors of Ref. Q 
carried out tri-dimensional Monte-Carlo simulations, in- 
cluding pion production effects and scattering off mag- 
netic field inhomogeneities. In their work, pion produc- 
tion is treated as a continuous energy loss process and its 
stochastic nature is neglected. The magnetic field config- 
uration is simplified to an assembly of coexistent bubbles 
of randomly aligned dipole fields, with a diameter given 
by the coherence length of the field. 

To properly calculate the desired likelihood, one has 
to consider both pion production and deflection on cos- 
mic magnetic flelds. Thus we devised a Monte Carlo 
code, which improves on the treatment in Ref. |31| on 
the points discussed in the following. 



Pion production is treated in the appropriate stochas- 
tic way, as described in Ref. |^o| ; namely, the occurrence, 
and the nature of the secondary and its energy are ran- 
domly drawn from the relevant probability distributions. 
Since the fractional energy loss of a proton in a pair pro- 
duction event is about the ratio of the electron to the 
proton mass, pair production can be included into the 
equations of motion as a continuous energy loss term, for 
which we used the expressions in Ref. . 

The LSMF are described as gaussian random fields, 
with zero ensemble average, and a power spectrum given 
by (B2(fc)) oc fc"« for k < 2n/lc, and (B^ik)) = other- 
wise. The cutoff, Ic, characterizes the coherence length of 
the field. The phase and direction of the field are drawn 
from uniform random distributions. Note that for the 
time delay calculation only the field component perpen- 
dicular to the line-of-sight contributes. This model can 
thus be characterized by three parameters: the power- 
law index ub, the coherence scale Ic, and the r.m.s. field 
strength, given by B^^^ = V/{27r)^ J d^kB^(k). 



The precise structure of the LSMF is unknown, but 
estimates can be made for different cases. Extra-galactic 
fields with cosmological origins have ub — and Ic — 
1 Mpc jssj. For halo fields to cause notable delays, Ic ^ 
1 kpc. Note that for these large coherence lengths the 
dynamical range probed by UHECRs is not very sensitive 

to HB- 

In addition to the LSMF defined above, there are local- 
ized fields associated with galaxies, clusters of galaxies, 
and the large scale galaxy distribution. The small filling 
factors of galaxies, fc ~ 10"'', and clusters, fc ~ 10~^, 
justify neglecting their effect. If galactic and cluster out- 
flows are efficient polluters of the extra-galactic medium, 
their effect can be modeled by the extra-galactic field con- 
sidered above. The filling factor for such pollution fields 
was estimated to ~ 10% for a field strength considerably 
smaller than ^ fiG ||3^ . 

Fields associated with the galaxies large scale structure 
may significantly affect UHECR propagation. In fact, 
the arrival directions of two of the three pairs observed 
by AGASA (see Table are very close to the super- 
galactic plane; reinforcing the claim for a supergalactic 
plane correlation of UHECRs . As was subsequently 
pointed out in Ref. |^ , this correlation seems to be even 
stronger than expected if the origin of UHECRs were 
associated with the large-scale galaxy structure, as the 
supergalactic plane alone is not a good approximation 
of this latter. As a solution, it has been suggested, in 
Refs. 38|, that the possible existence of strong fields, 
with strengths up to /xG and coherence lengths in the 
Mpc range, aligned along the large-scale structure, could 
produce a focusing effect of UHECRs along the sheets 
and filaments of galaxies. However, since UHECR prop- 
agation would take place mainly along the strong coher- 
ent field component, effects on the time delay would be 
much smaller than those induced by an incoherent field 
of comparable strength. Within a first approximation, 
therefore, we can neglect such coherent strong field com- 
ponents in our present analysis, which rests primarily on 
the time delay effect. In other words, our analysis is sen- 
sitive mainly to the weaker random field components that 
are not aligned with the large-scale structure sheets and 
filaments. 

We also neglect the influence of the magnetic field in 
the disk of our Galaxy. Denoting its typical strength 
by Bd and its scale by 1^, the deflection angle 9e for 
a proton at energy E arriving from above or below the 
galactic disk is 



5.5" 



E Y' f h 
lOEeVy VI kpc 



Ba 



10-6 G 



sm( 



(1) 



where 9 is the angle between the field polarization and 
the arrival direction of the particle. The corresponding 
time delay, te — 0^1^/2, is small compared to the mea- 
surement period of a few years and can be neglected 
for _E ^ 10 EeV in our analysis, if Ig < 100 pc and 
Sd ^ 10^6 G, which are typical values for the disk field. 
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FIG. 1. Contour plot for a distribution of events originat- 
ing from a bursting source, as projected in the time-energy 
plane. The corresponding parameters are D = 30Mpc, 
Brms ~ 10^^^ G, n_B — 0, Ic — 1 Mpc; the differential index of 
the energy injection spectrum is —2.0. Deviations from the 
mean correlation te oc , for E > 70EeV, reflect pion pro- 
duction on the CMB. In order to show the statistics at high 
energies, 40 contours with logarithmic decrements of 0.15 in 
the logarithm to base 10, are shown. 



We also remark that none of the AGASA pair arrival di- 
rections lie close to the galactic plane (see Table ||). We 
will assume that the actual time delay is dominated either 
by the extra-galactic magnetic field or by an extended 
galactic halo field. Under this assumption, observed time 
dispersions contain information only about the extra- 
galactic magnetic field or the galactic halo field . 

Once the configuration of the LSMF is set, via Fast 
Fourier Transform methods, on a lattice with inter-cell 
separation lc/2 and typical size 64*^ or 128^, a sample 
of nucleons is propagated over a distance D in a ran- 
domly chosen direction. The integration is achieved via 
the usual fourth order adaptive step size Runge-Kutta 
integration (see, e.g., Ref. Q). At each time step, the 
nucleons are subjected to deflection by the local magnetic 
held, are tested against pion production and /3 decay if 
the nucleon is a neutron, and pair production energy loss 
is taken into account. 

Angle-time-energy images of a discrete source at a fixed 
location depend on the specific LSMF structure which 
the nucleons encounter during their propagation |]39|] . In 
principle, one has to simulate a discrete source by emit- 
ting particles isotropically from the origin, and keep- 
ing those particles that arrive on the sphere of radius 
D within a solid angle segment small compared to the 
squared r.m.s. defiection angle. This is crucial in the 
limit of small deflection, D0e <C ^c; this limit is also the 
most probable for typical extra-galactic magnetic fields 
(see Eq. (||) below, see also Ref. In the opposite 



limit, D6e ^ Ic, where nucleons following different paths 
experience different magnetic field configurations, it is 
sufficient to send the particles isotropically, and average 
over all events independently of their arrival location. 
The effect of these two limits on the resulting angle-time- 
energy image of the source was pointed out in Ref. . 
In the limit of small defiection, the scatter around the 
mean correlations, te oc E^"^ and 6e oc E^^, is expected 
to be negligible below the photopion production thresh- 
old, with IS.E/E ~ 1%. In the opposite fimit, D9e > 1^ 
the stochastic deflections of different nucleons off differ- 
ent magnetic inhomogeneities imply a significant scatter, 
with /S.E/E ~ 30%. We confirmed that our code yields 
ditributions in energy and time that are consistent with 
the analytical expression given in Ref. ||2^] for this case, 
below the GZK cutoff. 

In order to simulate the small defiection limit, one 
has to restrain oneself to defining a receiving cone angle 
smaller than the r.m.s. defiection angle, and to selecting 
an emission cone angle not much larger than the receiv- 
ing angle. This will in general introduce a bias in the 
simulated energy spectrum, notably because the r.m.s. 
deflection angle is a function of energy. In order to cor- 
rect for this bias, we renormalize the time-integrated flux 
at a given energy to the corresponding flux for vanishing 
magnetic fields. This is appropriate in our case, where 
delay times are much smaller than energy loss distances, 
i.e., where the magnetic field has no infiuence on the time- 
integrated energy spectrum. Fig. |^ shows an example for 
the resulting distributions of arriving particles in time 
and energy for a bursting source, in the small defiection 
limit. For a more general discussion of the angle-time- 
energy images of UHECR sources, see Ref. [|9|. 



B. Likelihood 

We wish to evaluate the likelihood of the observations 
of the three AGASA pairs as a function of the different 
parameters characterizing the origin of UHECRs and the 
nature of their energy spectrum, as well as those charac- 
terizing the LSMF. We thus consider the following. We 
assume that the source emits UHECRs on the time scale 
Ts; for a burst, Ts ^ lyr. We group together the dis- 
tance to the source, Z?, the coherence length 1^ the r.m.s. 
magnetic field strength Brms, and the magnetic power 
spectrum index ns, in the average time delay te] this 
latter is calculated from the expressions given in Ref. , 
as: 



TE 



0.22 



riB 



D 



10-11 G 



30 Mpc 

Ir 



E 



100 EeV 



IMpc 



(2) 



This follows from a random walk argument, and is, 
strictly speaking, a good approximation only below the 
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pion production threshold and in the large deflection 
limit, DOe ^ Ic- We also consider the distance D as a 
separate parameter, as it governs the amplitude of pion 
production, hence the propagated UHECR energy spec- 
trum. 

For the total energy output of the source in UHECRs 
above an energy E, we assume a power law, £{> E) oc 
E~^~^'^, and consider three different representative val- 
ues for the index 7, namely 7 = 1.5, 2.0, and 2.5. A hard 
spectrum with 7 = 1.5 is typical for defect models |p3| , p5| , 
whereas softer spectra with 7 > 2.0 are typical for accel- 
eration at astrophysical shocks . With respect to ob- 
servations, it is more convenient to characterize the total 
energy output with the total number of particles N{Eo), 
that would be detected by the experiment with an energy 
E > Eq, over an infinite integration time. 

Here, /d denotes the experimental duty cycle with re- 
spect to having the source in the field of view, A is the 
effective detection area of the experiment under consider- 
ation, and fs{Eo) is a dimensionless factor of order unity 
that depends weakly on the energy injection spectrum 
cutoff. Eq. (^) is valid below the pion production thresh- 
old and in the following we take Eo — 10^'^EeV which 
is the energy above which the 3 AGASA pairs have been 
observed, within the energy resolution. We consider val- 
ues for iV(10^'^ EeV) between 1 and 4000 particles, which 
is a typical range to be expected for the fiducial values 
in Eq. (I). 

The main parameters defining the likelihood are there- 
fore Ts, rioo = TB^iooEcV, D, 7, and iVo = N{Eo = 
10^'^ EeV); secondary parameters are Ic and ns- The 
likelihood is calculated in the standard way for each ob- 
served event cluster denoted by zp, using Poisson statis- 
tics, 

where pj is the predicted number of events in cell j, and 
n(j, ip) is the number of observed events (either or 1 in 
our case) in cell j for the observed event cluster ip . Each 
cell is defined by a time coordinate and an energy. We 
binned the time-energy histogram to logarithmic energy 
bins of size 0.1 in the logarithm to base 10, and to 0.1 yr in 
linear time bins (we are limited to 0.1 yr time resolution 
because of memory size). The product in the formula 
above extends over all energy bins (from 10^'^ EeV to 
10^ EeV) and over all time bins that are comprised within 
an observational time window of length Tobs; we took 
Tobs — 5yr, as AGASA started operating in 1990, and 
the latest data reported date back to the end of October 
1995. 



The likelihood as given above has already been 
marginalized over two auxiliary parameters: a zero-point 
in time to, which defines the position of the observational 
time window on the time delay histogram of the UHE- 
CRs, and the path that the particle followed through 
the magnetic field. The first marginalization over to is 
carried out by drawing to at random a large number of 
times from a uniform distribution, calculating the likeli- 
hood, and averaging the resulting likelihood over to with 
equal weights. The second marginalization is carried out 
in a similar spirit, i.e., by separately calculating likeli- 
hoods for different realizations of the LSMF between the 
source and the observer and averaging them with equal 
weights. 

The likelihood defined in Eq. (^ is thus calculated in 
the following way: for each parameter (except Tg, 7, 
and No'- see below), and before each marginalization, 
typically 2 x 10^ particles are propagated. The result- 
ing histogram in time and energy is smeared out in en- 
ergy with the AGASA energy resolution AE/E ~ 30%; 
this histogram is also convolved on the time axis with 
a top-hat of width Ts to represent a uniform emission 
of particles during the time Ts. With the same popu- 
lation of propagated nucleons, one can construct differ- 
ent histograms depending on 7 and No- The nucleons 
are indeed injected with an artificially flat spectrum be- 
tween 10^'^ EeV and 10"* EeV, corresponding to 7 = 
to preserve satisfying statistics at high energies, and the 
counts on the detector are weighted in such a way as 
to reconstruct an injection spectrum with index 7. The 
likelihood is then marginalized over to, and calculated 
for each Tg. New samples of nucleons are propagated 
to marginalize the likelihood over the different magnetic 
field realizations; finally, the last loop closes, the above 
process is repeated for each value of rioo, and D. Due to 
the multi-dimensional nature of this likelihood, the size 
of the parameter space, notably the possible distances to 
the source, and the precision required (lyr corresponds 
to 3 X 10"^ Mpc) , these simulations are very time consum- 
ing. They were carried out on IBM at the Max-Planck 
Institut fiir Physik, Miinchen (Germany), and on Alpha 
at the Institut d'Astrophysique de Paris, Paris (France). 

We probe values for rioo between ~ 0.1 yr and 
1000 yr. The lower bound is given by the size of 
our time bins whereas the upper bound is motivated 
by Faraday rotation limits on the extra-galactic mag- 
netic field (e.g., Ref. |2^), which can be written as 
(B„ns/10-9G)(Ze/lMpcV/2 < 1, or, 

^-^^•^^^^^(so^)'^'- 

[see Eq. (|])]. A comparable bound comes from the re- 
quirement that the relative deflection of the events in the 
AGASA pairs was smaller than the angular resolution 
^1.6° 11. 

As to the source activity time scale, we choose the 
range 0.01 yr < Ts < lO'^ yr which, together with 1 < 
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Nq < 4000, brackets the range where the hkehhood func- 
tion should peak, given a couple of events observed within 
Syr. We note that, as long as the total UHECR output 
of the source, characterized by £{> 10^'^EeV) or Nq [see 
Eq. (^], is unconstrained, the likelihood analysis cannot 
distinguish between different values for Tg and tiqo that 
are considerably larger than the integration time Tobs 
of the experiment. In this limit, the likelihood becomes 
degenerate and only depends on the ratio(s) Nq/Ts or 
Nq/tiqq, i.e., the actual fluxes as seen by the experiment. 
Strictly speaking, it would be more appropriate to intro- 
duce the deviation around rioo in the latter ratio, instead 
of Tioo; however, at a given energy, this deviation is al- 
ways a finite fraction of rioo whose size only depends on 
whether the deflection is large or small compared to the 
coherence length [ p6| . 

It is interesting to note that pairs 2 and 3 do not show 
arrival energies inversely correlated to their arrival times, 
as would be expected from Eq. (||), whereas pair 1 does. 
As mentioned in the introduction, this flip around of time 
and energy for these pairs could be explained by the finite 
width of the te — E correlation, i.e., the combination of 
the scatter due to the stochastic nature of the magnetic 
field, and of the finite instrumental energy resolution. 
Nevertheless, this could still be insufficient for pairs 2 
and 3; in this case, these pairs would tend to favor a 
continuous source, where a flip around can simply result 
from different emission times. We will observe each of 
these tendencies in Sec. 4. 

Finally, we restrict ourselves to the limit of small de- 
flections, D9e ^ ^c- In effect, since the three pairs ob- 
served by AGASA imply Be < 2.5° ~ 0.044, the limit 
D9e ^ would require D/lc ^ 25 and thus a rather 
large lattice in our code. Moreover, since the AGASA 
pairs also imply D 50 Mpc, due to pion production, 
this case would only be realized if Ic ^ 2 Mpc. This 
is not likely in the case of an extra-galactic magnetic 
fleld, due to possible damping of the fleld on small scales 
< 1 Mpc [||. If the deflection of UHECRs is dominated 
by the halo magnetic fleld, the distance to the source 
in Eq. has to be replaced by the linear extension of 
the magnetic halo Zh (if the source lies outside of this 
halo); in this case, the limit D9e ^ ^c, would necessitate 
Ic < 4kpc if In < 100 kpc. 



III. RESULTS 

Since the likelihood function depends on so many vari- 
ables we will present several one- and two-dimensional 
cuts through the parameter space. Because the values 
of most of the parameters are basically unconstrained it 
is also convenient to marginalize over one or several of 
the variables. In order to do so, one has to choose a 
suitable Bayesian prior. We will always assume a Jeffrey 
prior, i.e., uniform distributions in the logarithm of the 
variables tiqo, Is, and A^o- This choice is motivated by 
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FIG. 2. Contour plots of the log-likelihood logj^o i^i 
Tioo — Ts plane for pair 1. The corresponding parameters are 
D — 30 Mpc, ub = 0, Ic = 1 Mpc; A'o and 7 as indicated. The 
bottom panel corresponds to the maximum of the likelihood. 
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FIG. 3. Same as Fig. 2, but for pair 2. 
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FIG. 4. Same as Fig. 2, but for pair 3. 

the "scale invariance" of the problem, namely the fact 
that the time delay depends on powers of the magnetic 
field strength and coherence length, which are totally un- 
known, as well as the fact that the emission time scale 
and energy output are unknown within many orders of 
magnitude. We will thus consider the following types of 
representations in all of which D — 30Mpc, rts = 0, and 
Ic — IMpc take specific constant values, for illustration: 

1. For a given Nq, contour plots of the likelihood 
\ogiQ{C) are shown in the tiqo— Is plane, for two ex- 
treme values of the power-law index of the injection 
spectrum, 7 = 1.5, and 7 = 2.5, see Figs. (||— §)■ 

2. The likelihood C{No), as marginalized over Ts and 
Tioo, is shown in Fig. (||). Via Eq. (||), this de- 
termines the most likely value for the total energy 
output £{> lOEeV) of the source. 

3. The likeHhood C{tioo,Ts <^ O.lyr) marginalized 
over A'o is shown in Fig. (||). This particular 
marginalization assumes a bursting source, and 
thus allows to determine the most probable corre- 
sponding time delay. Notably, the model of origin 
of UHECRs in cosmological 7-ray bursts requires 
rioo>50yr 0,0. 

4. The likelihood jC{Ts), as marginalized over tiqo and 
A^'o, is shown in Fig. (Q). This marginalization al- 
lows to test the respective significance of a bursting 
or a continuous source. 

From the figures shown here it is obvious that there 
are too many parameters for the current number of ob- 
servations available. However, the following trends are 
discernible: 



• In the limit of a high number of particles, the typ- 
ical time dispersion, i.e., either Ts or the deviation 
around tiqo, scales like No, so that the particles are 
spread out in time on the detector; i.e., so that the 
observed flux matches the injected flux. The like- 
lihood thus becomes degenerate, and depends only 
on the ratios Nq/Ts, and/or Nq/tioq, as argued in 
Sec. 2. This degeneracy arises whenever Ts ^ TLbs, 
or the dispersion of the time delay for the lower en- 
ergy particle in the pair, is much larger than Tobs- 
This can be seen best in Figs. 18[ In fact, pairs 2 
and 3 (see Table |), in which the higher energy 
particle arrived later, favor, respectively marginally 
and clearly, a comparatively large Tg. 

• Pair 1 marginally (i.e., within factors of a few in 
the likelihood) favors a hard injection spectrum, 
which is expected in order to produce a 200 EeV 
event. For tiqo ^ a few years it allows for a burst- 
ing source. This is mainly due to the fact that its 
energies and arrival times are inversely correlated, 
as expected from Eq. (p|). The small magnitude of 
the corresponding tiqo S 1 yr results from the small 
difference in arrival times as compared to the signif- 
icant difference in arrival energies (see Fig. ^). The 
origin of this pair is therefore not straightforward 
to account for. 

• Pair 2 is insensitive to 7, since the energies of the 
particles are very close to each other. It allows for 
a burst with a large time delay, although the arrival 
times and energies are not inversely correlated. As 
argued previously, the flip around of the energies 
can be explained by the finite instrumental energy 
resolution. In our case, the instrumental resolution 
AE/E ~ 30% dominates over the intrinsic scat- 
ter due to the scattering off the magnetic field. 
This width introduces a dispersion in time delay, 
^te/te = 2AE/E. Hence, provided the time de- 
lay is sufficiently large, there is significant scatter 
both in energies and arrival times. This explains 
the tendency of pair 2 to accept a large time delay. 
Although the marginalized likelihood C{Ts) slightly 
favors Tg 10 yr, this is not very significant and 
still consistent with a burst (see Fig. |^). This pair is 
thus consistent with an origin in cosmological 7-ray 
bursts. 

• Pair 3, however, appears to be completely inconsis- 
tent with a burst for time delays rioo ^ 100 yr. For 
larger time delays and correspondingly high total 
expected number of particles Nq, this pair could 
be consistent with a burst. However, such high 
time delays might in turn be inconsistent with the 
Faraday rotation limit, as given by Eq. (^). Fur- 
thermore, the likelihood is suppressed by roughly 
a factor of 100 for such a scenario (see Figs. ^ 
and 1^). Overall, and within the parameter space 
probed here, this pair seems to favor an origin in a 
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No 

FIG. 5. Likelihood marginalized with respect to rioo and 
Ts, plotted vs No, for pair 1 (bottom panel), pair 2 (middle 
panel), and pair 3 (top panel). The solid lines correspond to 
7 = 1.5, the dotted lines to 7 = 2.0, and the dashed lines 
to 7 = 2.5. The parameters D, hb, and Ic are the same as 
in Figs. 2—4. The decrease for A^o > 500 for pair 1 and 3 is 
explained by the fact that in this range one would expect more 
than 2 events per 5 yr since Ts < 10^ yr in our simulations. 
Due to the lower energies this can be compensated by a large 
time delay In case of pair 2. 

continuous source with a high total energy output, 
corresponding to Nq > 100. 

• In all cases the likelihood function tends to peak in 
the range where rioo ^ 10 yr. This is not significant 
for pair 2, which seems consistent with a burst and 
a comparatively long time delay; for pairs 1 and 3, 
however, the likelihood peak is at least a factor 10 
higher than any value in the domain tiqo > 10 yr. 
Using Eq. (|), this can be transformed into the ten- 
tative limit 

Sons < 6 X 10-" ( ( — ' G . 

^ \lMpcJ \30MpcJ 

(6) 

If confirmed, or possibly even improved by future 
data, this constraint would be more stringent than 
existing limits [ p7| by more than an order of mag- 
nitude. 

We also performed runs for D = lOMpc and D = 
60Mpc, as well as for ub = 4 and for 1^ = lOMpc. For 
corresponding values of the parameters in Eq. (^, the 
likelihood does not significantly depend on and Ic- 
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FIG. 6. Likelihood marginalized with respect to A'o, as- 
suming a bursting source, Ts <^ l-yr, plotted vs rioo, for 
pair 1 (bottom panel), pair 2 (middle panel), and pair 3 (top 
panel). The solid lines correspond to 7 = 1.5, the dotted lines 
to 7 = 2.0, and the dashed lines to 7 = 2.5. The parameters 
D, ns, and /c are the same as in Figs. 2—4. This particular 
marginalization allows to estimate the most probable time 
delay, assuming a bursting source. 
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FIG. 7. Likelihood marginalized with respect to rioo and 
TVo, plotted vs Ts, for pair 1 (bottom panel), pair 2 (middle 
panel), and pair 3 (top panel). The solid lines correspond to 
7 — 1.5, the dotted lines to 7 = 2.0, and the dashed lines to 
7 — 2.5. The parameters D, ub, and Ic are the same as in 
Figs. 2-4. 
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excluded for the sources of the three AGASA pairs. 




0.1 1,0 10.0 100.0 1000.0 

Tg [yr] 

FIG. 8. Same as Fig. 7, but for a source at D — 100 kpc and 
a galactic halo field with ns = 4, and Ic — 30 kpc, extending 
over 100 kpc. 

For D ^ 60 Mpc the dependence on D is also insignifi- 
cant, except possibly for pair 1 which favors D < 30 Mpc 
within a factor 3 — 4 in the likelihood. 

Furthermore, we performed simulations for sources lo- 
cated within an extended galactic halo, in which case 
the time delay is dominated by the halo magnetic field. 
We chose a halo extension of 100 kpc with a coherence 
scale Ic = 30 kpc, = 4, although the results again 
are not sensitive to these parameters. The maxima of 
the likelihood do not significantly differ from those ob- 
tained previously for a dominant extra-galactic magnetic 
field. The main difference between these two limits is 
the predominance of a soft spectrum 7 ~ 2.0 for pair 1, 
in the halo field case, even though this pair includes the 
200 EeV event. This is demonstrated in Fig. || and can 
be attributed to the absence of frequent pion production. 
Furthermore, if D is now interpreted as the minimum of 
the source distance and the extension of the magnetized 
halo, Eq. (||) translates into the tentative constraint on 
the halo field 

Finally, we mention that the results presented in this 
section do not depend significantly on the lower bound- 
ary Eq of the energy range over which we compute the 
likelihood function Eq. (Q) as long as it is not much 
lower than the lowest event energy in the pairs, within 
the energy resolution. Lowering Eq tends to make the 
7— dependence of the likelihood function more significant, 
but neither a hard nor a soft injection spectrum can be 



IV. DISCUSSION AND CONCLUSIONS 

Under our assumptions on the configuration of the 
LSMF, namely a power-law power spectrum down to 
some cutoff scale Mpc and ^ 1 kpc for time de- 

lays dominated by the extra-galactic or the galactic halo 
field, respectively, and within the limited range probed 
for the average time delay tiqo , our analysis suggests that 
all three AGASA pairs originated in different sources! 
Pair 1 would originate from a burst with a short time de- 
lay, Tioo lyr, i.e., a time delay inconsistent with that 
required in a cosmological GRB model. Pair 2 could orig- 
inate from a burst, if the average time delay is sufficiently 
large, tioo ^ 50yr, as advocated in the cosmological GRB 
model, although it marginally favors a continuously emit- 
ting source with a high energy output. Finally, pair 3 
cannot originate from a burst, unless the time delay is 
extremely large, and probably too large when compared 
to the Faraday rotation bound; it has to originate within 
a continuously emitting source, with Tg > lOOyr, and 
A^o ^ 100, in agreement with conventional sources such 
as powerful radio-galaxies, for instance. 

Let us now address the possibility of bursting topo- 
logical defect sources. For example, certain classes of 
cosmic string loops might collapse and release all of their 
energy in form of UHECRs within about one light cross- 
ing time Tg <^ 1 yr |42j , for symmetry breaking scales 
near the grand unified scale and typical total burst ener- 
gies ~ 10^^ erg. For an LSMF > 10~^^G, events above 
~ 80 EeV are predicted to be most likely 7-rays, whereas 
around 50 EeV an approximately equal amount of pro- 
tons is expected To take into account the likely 
7— ray nature of the higher energy events in these models, 
a more accurate treatment should include propagation 
and deflection of electromagnetic cascades in the large- 
scale magnetic field. However, since for a given energy 
the amount by which an electromagnetic cascade parti- 
cle and a nucleon is deflected and delayed is comparable, 
the topological defect scenario for the origin of UHECRs 
is not inconsistent with our above results, notably as far 
as pair 1 is concerned. In this respect, we note that the 
muon contents of the AGASA pairs air showers are not in 
contradiction with interpreting the higher energy event 
as a 7-ray j^] . The most efficient test of the topological 
defect scenario clearly resides in the chemical composi- 
tion at energies > 80EeV, and improved data on UHECR 
composition could rule out this model in the future. 

We do not wish to emphasize too strongly the above 
results, as quite a number of approximations had to be 
made, that restrict the range of parameter space available 
to us, and, furthermore, the statistics of the AGASA ob- 
servations are limited. The main objective of the present 
work was, rather, to introduce the numerical code devel- 
oped, and to give an example of what could be achieved 
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by future detectors. From our experience, it seems rea- 
sonable to say that one or several clusters of more than 

10 particles per source, would allow to pin down most 
of the parameters, notably: the magnetic field strength 
-Brms and coherence length Ic, the distance to the source 
D, the emission time scale Tg, and the time delay tioo- 
Other parameters such as the differential energy index 7 
and the index of the power spectrum of magnetic inho- 
mogeneities ns (if any), are more difficult to grasp. In 
a separate paper, we present qualitative aspects of the 
angle-time-energy images of UHECR sources the 
qualitative features given there would allow to consid- 
erably restrain the size of the parameter space, and even 
give estimates of some of the above parameters. This 
would constitute a very useful first approach to a full 
likelihood evaluation. 

In Ref. Q] it was shown that the 7— ray flux above 
10^^ eV contains information about the strength of the 
extra-galactic magnetic field. Here, we have shown 
that a likelihood analysis of clusters of charged UHE- 
CRs from a common extra-galactic source is also sensi- 
tive to the strength and the coherence scale in a range 
lO-i^G < B,^s ^ lO-^G, with > IMpc, and in 
addition to a possible galactic halo field in the range 
10-« G < B,r,,s < 10"^ G, with k > 1 kpc. The deflection 
and delay of such UHECR clusters can thus be used as a 
new tool to probe the LSMF in a regime that is interme- 
diate between field strengths accessible to conventional 
methods such as measuring the Faraday rotation of po- 
larized light (see, e.g., Ref. |^), and much weaker fields 
that could be detected by observing electromagnetic cas- 
cades in the TeV range The UHECR pairs ob- 
served by AGASA already indicate a trend towards com- 
paratively small time delays, that can be translated into 
a tentative upper limit on the magnetic field strength, 
roughly an order of magnitude better than the Faraday 
rotation limit. 

Future instruments in construction or in the proposal 
stage such as the Japanese Telescope Array the 
High Resolution Fly's Eye |^^, and the Pierre Auger 
Project 1^ will have the potential to test whether there 
is significant clustering of UHECRs. The latter experi- 
ment, with an angular resolution of a fraction of 1° and 
an energy resolution of ~ 10%, should detect clusters of 
20 — 50 events if the clustering observed by AGASA is 
real. With such statistics, it will be possible to answer 
the question on the nature of the sources of such UHECR 
clusters, and, in turn, to use these as "candles" to probe 
cosmic magnetic fields. 
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